function convInd = checkConvex(restr,phi,opt,convInfo)
% Check whether identified set (IS) is convex using Proposition 3.
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options
% - confInf: structure containing information about convexity of IS

F = restr.F;
f = restr.f;
S = restr.S;
s = restr.s;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
jshock = opt.jshock;
signNorm = opt.signNorm;                         

N = size(Sigmatr,1); % No. of variables in VAR

% convInd = 0 if Proposition 3 cannot guarantee convexity of IS.
convInd = 0; 

if convInfo.propCase == 3 || convInfo.propCase == 5
        
    % Use Proposition 3(I)(iii). That is, check if there exists
    % 1 <= istar <= (jshock-1) such that [q_1,...,q_istar] is exactly 
    % identified and f_i < n-i for all i = istar+1,...,jshock.
          
    % First consider istar = 1. q_1 can only be exactly identified if
    % f_1 = N-1 and rank(F_1(phi)) = f(1). For 3(I)(iii) to be 
    % satisfied, also require  f_i < N - i for all i = 2,...,jshock.
        
    if (f(1) == N-1) && (rank(F(1:f(1),:)) == f(1)) && ...
            (all(f(2:jshock) < N - (2:jshock)'))

        if convInfo.propCase == 3
            
            % istar = 1 and identified set is convex.
            convInd = 1;
            return
            
        elseif convInfo.propCase == 5
            
            % Add check for vector satisfying sign restrictions!
            
        end

    end
        
    for istar = 2:(jshock-1) % Consider 2 <= istar <= jshock - 1

        % [q_1,...,q_istar] can only be exactly identified if 
        % f_i = N-i and rank(F_i(phi)) = f(i) for i = 1,...,istar.
        % For 3(I)(iii) to be satisfied, also require  
        % f_i < N - i for all i = istar+1,...,jshock.

        rankF = zeros(istar,1);
        rankF(1) = rank(F(1:f(1),:));

        for jj = 2:istar

            rankF(jj) = rank(F(sum(f(1:jj-1))+1:sum(f(1:jj)),:));

        end

        if all((rankF == f(1:istar))) && ...
                all(f(1:istar) == (N-(1:istar)')) && ...
                all(f(istar+1:jshock) < (N-(istar+1:jshock)'))

            % Assess whether Algorithm 3 yields a unique set of
            % orthonormal vectors [q_1,...,q_istar].
            rankFq = zeros(istar,1);
            signNormCheck = zeros(istar,1);

            F1 = F(1:f(1),:);
            qtilde = null(F1); % q_1

            if signNorm == 0 % Sign normalisation on A0

                signNormCheck(1) = Sigmatrinv(:,1)'*qtilde;

            else % Sign normalisation on A0^(-1)

                signNormCheck(1) = Sigmatr(1,:)*qtilde;

            end

            for jj = 2:istar

                Fi = F(sum(f(1:jj-1))+1:sum(f(1:jj)),:);
                rankFq(jj) = rank([Fi; qtilde']);
                q_i = null([Fi; qtilde']); 

                if signNorm == 0 % Sign normalisation on A0

                    signNormCheck(jj) = Sigmatrinv(:,jj)'*q_i;

                else % Sign normalisation on A0^(-1)

                    signNormCheck(jj) = Sigmatr(jj,:)*q_i;

                end       

                qtilde = [qtilde; q_i];

            end

            if all(rankFq(2:istar) == (N-1)) && ...
                            all(signNormCheck~=0)
            
              % [q_1,...,q_istar] are exactly identified.
              
              if convInfo.propCase == 3
              
                  convInd = 1; % Identified set is convex by 3(I)(iii).
                  return
                  
              elseif convInfo.propCase == 5
                  
                  % Use 3(II)(v).
                  Fjstar = F(sum(f(1:jshock-1))+1:sum(f(1:jshock)),:);
                  qstar = null([Fjstar; qtilde]);
                  
                  % Construct matrix of inequality restrictions on q_jshock
                  % (including sign normalisation).
                  Sjstar = zeros(s(jshock)+1,N);
                  Sjstar(1:s(jshock),:) = ...
                      S(sum(s(1:jshock-1))+1:sum(s(1:jshock)),:);

                  if signNorm == 0 % Sign normalisation on A0            

                      Sjstar(end,:) = Sigmatrinv(:,jshock)';

                  else % Sign normalisation on A0^(-1)

                      Sjstar(end,:) = Sigmatr(jshock,:);

                  end

                  for kk = 1:size(qstar,2)  % May be multiple possible vectors

                      if all(Sjstar*qstar(:,kk) > 0) || ...
                              all(-Sjstar*qstar(:,kk) > 0)

                          convInd = 1; % Identified set is convex
                          return

                      end

                  end
  
              end

            end

        end     
        
    end

elseif convInfo.propCase == 4
    
    % Use Proposition 3(II)(iv) to check convexity of IS.

    % Construct matrix of equality restrictions on q_jshock.
    Fjstar = F(sum(f(1:jshock-1))+1:sum(f(1:jshock)),:);

    % Construct matrix of inequality restrictions on q_jshock (including
    % sign normalisation).
    Sjstar = zeros(s(jshock)+1,N);
    Sjstar(1:s(jshock),:) = S(sum(s(1:jshock-1))+1:sum(s(1:jshock)),:);
    
    if signNorm == 0 % Sign normalisation on A0            

        Sjstar(end,:) = Sigmatrinv(:,jshock)';

    else % Sign normalisation on A0^(-1)

        Sjstar(end,:) = Sigmatr(jshock,:);

    end

    % Check whether there exists a unit length vector satisfying 
    % the equality restrictions in Fjstar and the sign 
    % restrictions (and normalisation) in Sjstar.
    % First compute orthonormal basis for nullspace of Fjstar (i.e. a 
    % unit length vector satisfying equality restrictions in Fjstar).
    qjstar = null(Fjstar); 

    for jj = 1:size(qjstar,2)  % May be multiple possible vectors

        if all(Sjstar*qjstar(:,jj) > 0) || all(-Sjstar*qjstar(:,jj) > 0)

            convInd = 1; % Identified set is convex

        end

    end

end

end